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ABSTRACT 

A three-dimensional analytical model for large water bodies is presented. 
Time histories and spatial distribution of pressure, velocity and tempera- 
ture in water bodies, subjected to thermal discharge, are determined em- 
ploying a digital computer. The dynamic response is obtained for a 
rectangular water body by applying a finite difference method to the 
mass, momentum, and energy balance equations. The distinctive feature of 
this analysis, as compared to previous studies, is the calculation of 
pressure and water level from equations of motion without simplifying 
assumptions such as hydrostatic pressure approximation and rigid-lid con- 
cept. The mathematical formulation is verified by applying this analysis 
to cases where the final steady state flow patterns have been determined 
analytically or experimentally by others. The problem of thermal dis- 
charge entering a river is then analyzed. The time histories of the 
velocity and temperature distribution are obtained. These results pro- 
vide the values of temperature rise and the rate of temperature rise 
needed for the assessment of thermal pollution in water bodies. 


INTRODUCTION 

Thermal power and industrial processing plants use large quantities of 
cooling water from natural or artificial water bodies and return the 
water to the source at higher temperatures. When the thermal discharge 
from a plant, into a water body, raises the water temperature by such an 
extent as to damage aquatic life or other legitimate uses of the water 
source, some degree of thermal pollution exists. 

To assess the extent of a thermal pollution, utilities, ecologists, fed- 
eral and state agencies as well as the public are Interested in determin- 
ing the space and time distribution of temperature within a water body, 
when it is subjected to a prescribed thermal load. 

The theoretical analysis of the dynamic behavior of a large water body, 
subjected to thermal discharge, is based on the solution of space- and 
time-dependent conservation and state equations. Conservation of mass, 
momentum, and energy for the water flow in the Impoundment together with 
the equation of state, provide a sufficient number of partial differen- 
tial equations and algebraic relations in space and time for the solution 
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of the problem. Although, these equations have been known for over a 
century, due to their highly nonlinear nature, a direct closed form 
analytical solution is considered practically impossible for most general 
cases involving realistic geometries and boundary conditions. The avail- 
able mathematical solutions of problems involving thermal discharge ate 
severely limited by the many simplifying assumptions made in order to 
achieve a solution. In fact, analytical solutions to the problems of 
thermal discharge are obtained only for problems involving simple geom- 
etry and boundary conditions. Solutions to more realistic geometries and 
boundary conditions can be sought only by the development of computer 
modeling in three-dimensions under transient conditions. 

Current three-dimensional computer simulation of thermal discharge into 
water bodies are based on a number of simplifying assumptions discussed 
hereunder: 

1. Compliance with the Conservation of Mass and Momentum 

Generally, the three components of water velocity are found from the lon- 
gitudinal, lateral, and vertical momentum equations. These three compon- 
ents of velocity should satisfy the conservation of mass equation. A 
major difficulty arises when simultaneous compliance with the principles 
of conservation of mass and momentum cannot be ensured. Under the 
Boussinesq approximation, the conservation of mass reduces to zero div- 
ergence for the water velocity (with no mass storage term) i.e. the sum 
of all inflows to and outflows from any fluid control volume in the flow 
field must vanish. Since the new velocity components, calculated by the 
integrations of the momentum equations do not necessarily satisfy the 
continuity equation, the mass balance would be disturbed, thereby yield- 
ing a small surplus or deficit inflow at certain control volumes. To 
resolve this dilemma, Brady and Geyer [1] hypothesize that the surplus or 
deficit inflows redistribute in all directions by equal magnitude. Ob- 
viously, the disadvantage of this velocity adjustment technique is that 
the adjusted velocities no longer satisfy the momentum equations. Con- 
sidering a water body, such as a river, flowing longitudinally with a 
laterally uniform velocity distribution, and applying the appropriate 
boundary conditions on the lateral velocities at the thermal discharge 
location, they found that the out-of-balance surplus inflow to each ad- 
jacent cell is exceptionally large. When their model attempts to distri- 
bute this surplus inflow among the surrounding fluid cells, the resulting 
vertical velocity adjustment becomes excessively large which aborts the 
run. Brady and Geyer, referred to this problem as the "vertical over- 
responsiveness of the model" and devised schemes to suppress the associ- 
ated undesirable vertical fluctuations by either: 1) temporarily sup- 

pressing vertical redistribution such that each fluid layer achieves its 
own mass balance independently; or 2) temporarily enforcing a rigid lid 
at the fluid free surface such that a zero vertical velocity is main- 
tained at the most upper layer of the model at every fluid cell. The 
rigid lid assumption is used by many authors to simplify the solution to 
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the stratification and circulation in water bodies [2, 3, 4]. Obviously, 
this situation calls for an improved approach for a simultaneous compli- 
ance with the conservation equations of mass and momentum. 

2. Surface Flow Distribution 


The treatment of the slope variation at the water free surface consti- 
tutes another major difficulty. The upwelling (or downwelling) flow at 
the surface creates positive (or negative) surge waves [5, 6]. The mathe 
matical modeling of these phenomena under unsteady flow conditions is 
complicated. To resolve this difficulty Waldrop and Farmer [7] first as- 
sumed that the vertical component of velocity near the surface will be 
redirected in the horizontal direction (in x and y directions) and the 
surface will move. Later Waldrop and Farmer [8] in a major effort to re- 
solve the surface flow distribution problem introduced a mass balance at 
the free surface and extended their horizontal momentum equations to ele- 
ments near the water surface. This situation calls for further studies 
of the behavior of the free surface under unsteady flow conditions. 

3 . Pressure Distrxbution 

The velocity components in the flow field are extremely sensitive to 
values of nodal pressures. Slight changes in the pressure distribution 
would affect the circulation in the water body considerably. The calcul- 
ation of the correct pressure distribution is, therefore, of paramount 
importance. Most authors [9, 10, 11, 12, 13, 14, 15] assume that the 
total dynamic pressure at each point in the flow is equal to hydrostatic 
pressure obtained under static conditions. This assumption, referred to 
as hydrostatic approximation, leads to inaccuracies in regions of severe 
upwelling and downwelling which will adversely affect the correctness of 
circulation velocities. Thus, a better technique for the computation of 
the pressure field is called for. 

4. Numerical Stability 

Calculations of space and time increments for obtaining a stable numeri- 
cal solution of partial differential equations of mass, momentum and 
energy conservation is a difficult task. Most authors resort to overly 
simplified criteria for the establishment of the upper bound of the inte- 
gration time step for an assumed set of space increments . Brady and 
Geyer [1] state that the maximum integration step appears to be limited 
by relationships between the velocities and the smallest cell dimensions 
in each direction. Waldrop and Farmer [7], and Harlow and Welch [16] re- 
sort to a stability limit calculated on the basis of one-dimensional, un- 
steady incompressible flow equations with a free surface. This limit 
imposes an upper bound on their integration time step equal to the ratio 
of the smallest cell size and the surface wave speed. More accurate 
means for the prediction of the space and time increments are needed. 

Most authors believe that much work is still required to refine various 
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features of the present unsteady three-dimensional computer models before 
the present models can be applied to situations involving significant 
stratification effects, or fluctuations in flow in water bodies. 

The main objectives of the present study are as follows: 

1) To develop an improved three-dimensional analytical model for the 
mathematical description of a large rectangular water body subjected to a 
thermal discharge. 

2) To develop a numerical integration technique for solving the above 
mathematical model on a digital computer. 

3) To develop a digital computer program for predicting the thermal 
stratification and circulation phenomena in a given rectangular water body 
subjected to thermal discharge and to determine the time histories and 
spatial distribution of pressure, velocity and temperature fields within 
the water body. 

The major contribution and distinctive features of the present study, as 
compared with the previous works in this general area, can be summarized 
as follows: 

A) Compliance with the conservation of mass and momentum is obtained by 
the introduction of two flow regions in the entire flow field: a) the 
water-level region containing a portion of water near the free surface in 
which the water level rises or falls, as the case may be, during the dy- 
namic solution; and b) the sub-water-level region located under the water 
level region which remains totally filled with fluid at all times during 
the transient. As will be seen later, the superimposition of a three- 
dimensional grid on these two regions facilitates the simultaneous com- 
pliance with the principles of conservation of mass and momentum for the 
entire flow field without any mass imbalance. 

B) A different set of differential equations for the conservation of 
mass , momentum and energy is applied to the cells located in the water- 
level and sub-water-level flow regions. This is necessary because the 
cells in the sub-water- level region have a variable height while the cells 
in the sub-water-level region have a constant height. 

cj The pressure distribution is obtained by combining the momentum and 
continuity equations. This feature eliminates the need for the hydro- 
static pressure approximation. 

D) A detailed numerical stability analysis is performed which provides 
accurate criteria for the selection of the space and time increments. 


MATHEMATICAL FORMULATION 


ANN 


-4- 



VI-B-41 


Simplifying Assumptions 

The mathematical formulation in this study is based on the following as- 
sumptions : 

1) The equations of motion, energy and continuity are applied in their 
time-smoothed form to turbulent incompressible three-dimensional flow with 
Cartesian coordinates, as shywn in Figure 1, with z positive upward. 

2) It is assumed that Boussinesq's approximation applies. In other 
words, densities are treated as constants except in terms involving grav- 
ity. This allows natural circulation to take place. Furthermore, den- 
sity in the body force term is considered to be a function of temperature 
only. 

3) The effects of turbulence are modeled by using eddy transport coeffi- 
cients . Horizontal and vertical momentum eddy viscosities and thermal 
eddy dif fusivities are considered as constants throughout the body of 
water though different magnitudes for horizontal and vertical directions. 

4) It is assumed that there are no internal heat sources and the heat ex- 
change between the water body and the atmosphere takes place near the 
water free surface. 

5) Loss of mass due to evporation at the surface and conductive heat 
transfer through the impoundment solid boundaries are generally small and 
are neglected. 

6) The Coriolis forces acting on the water body is considered to be neg- 
ligible. 

The boundary conditions of the problem involve; 1) the geometry of the 
Impoundment Including the thermal discharge inlet, as well as the river 
Inflows and outflows configurations; 2) the mass flow rates of the in- 
flows across the boundaries; and 3) the level, pressure and temperature 
along the inlet boundaries. It is assumed that all of the above para- 
meters are known. 

The initial conditions of the problem includes the pressure, velocity, 
level and temperature distributions of the impoundment at time t = 0. 

These variables are computed by the present program by first performing 
a dynamic analysis under no thermal discharge conditions and then using 
the calculated pressure, velocity, level and temperature distributions, 
as initial conditions, in a dynamic analysis involving a thermal dis- 
charge . 

Governing Equations 

The mathematical formulation of the sub-water-level flow region consists 
of the following fundamental equations [17]. 
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Momentum equations : 
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Equation of state: 
P = P (T) 


( 6 ) 


In this study, the thermal discharge parameters are used as reference 
values for the non-dimensionalization of the governing equations. Let 
U^, Pq, Tq and dQ be the velocity, density, temperature and half width of 

thermal discharge flow respectively. The dimensionless quantities are 
defined by 
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Employing the above dimensionless variables in the continuity, momentum 
and energy equations, one obtains the time derivatives for velocity com- 
ponents u, V, w and temperature T in parabolic differential form as 
follows: 
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In the above equations, the following dimensionless nximbers are used 
1) Froude number Fg = Ug/ (gdg)^^^ (17) 
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Reynolds number R = U-d„/v_ 
0 0 0 0 
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Equations (11) to (16) are a set of six equations that will be used for 
the determination of six unknowns, u*, v*, w*, p*, p* and T* as functions 
of time and space. 


The final equations, including the pressure equation, are obtained from 
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Mu I (juations as follows. First, the continuity equation is dif- 

torontiated with respect to time to give 
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The order of differentiation in equation (19) is then interchanged 
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( 20 ) 


To facilitate spatial differentiations indicated by eq. (20), eq. (11) 
through (13) are first rewritten in the following form 
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( 27 ) 


substituting eqs. 


2 * 

P = 


(21) through (23) into 
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This elliptic differential equation will provide a means for the calcula- 
tion of the pressure distribution in the water body. Eqs. (27), (24), 

(25), (26), (28), (21), (22), (14) and (16) constitute our final equa- 
tions which will be solved for updated values ofw,Q,Q,Q,p,u, 
v”, T*, and p* respectively. It should be noted that: " ^ 

The vertical momentum equation (23) is used indirectly in this computation 
through Its usage in the derivation of the pressure equation (28) . Since 
the water level variation provides a vertical freedom for the water body, 
It is imperative to calculate the vertical velocity component w from the 
continuity equation and not by the application of the momentum equation in 
the vertical direction equation (23). However, the momentum balance In 
the vertical direction is fully satisfied since eq. (23) is explicitly 
used in the derivation of the pressure equation (28). In fact, eq. (23) 
may be easily derived from the final equation stated above. 


The mathematical formulation for the water-level flow region is derived 
by the application of the fundamental equations of mass, momentum and 
energy as follows. The water level in cells located at the air-water 
interface varies with time and space. For this reason, the continuity 
equation, as given by eq. (27), is not valid for these cells. Applica- 
tion of a mass balance to such a cell with all sides fixed in space ex- 
cept the top which moves with the water level, yields 

(pAxAy H) = (puAy H)^ - (puAy H)^^^ (29) 

+ (pvAx H)y - (pvAx H)y^^y 
+ (pwAxAy)^ - (pwAxAy)^^g 


The last term represents the mass flux leaving the water level surface. 
Since the evaporation losses are considered insignificant, this term is 
equal to zero. Dividing both sides of the above equation by pAxAy, con- 
sidering P to be constant and equat to Pq consistent with Boussinesq ap- 
proximation used earlier, letting Ax, Ay approach zero and using 
dimensionless variables yields 
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where 


* H 

" -T. 


(31) 


Similarly, application of the energy and momentum balance to the water- 
level region gives: 
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The pressure at the air-water interface is always atmospheric and can be 
considered equal to zero no matter what the water level may be at any in- 
stant of time. For this reason, the pressure at the center of the water 
level can be considered to be hydrostatic, i.e. for the water level ele- 
ment 


it * * 

P = P H 


(38) 


This equation may also be derived by applying a momentum balance, in the 
vertical direction, in the water-level flow region and noting that the 
water vertical velocity in this region is relatively small and therefore 
negligible. The pressure distribution given by eq. (38) is substituted 
into the horizontal momentum equations. 


It Is important to realize the difference between the flow patterns in 
the flow field and the water level elements. In the flow field elements, 
an increase in the fluid boundary velocity, will induce horizontal velo- 
city components. However, the time constant associated with the momen- 
tum equations will not permit these equations to respond instantaneously 
to the boundary velocity disturbance at the river or thermal discharge 
inlets. Because of the flow incompressibility and since the fluid ver- 
tical inertia is considerably smaller than the horizontal inertias, the 
mass imbalance will result in a positive vertical fluid motion desig- 
nated as "upwelling" , This vertical motion will continue upward until it 
reaches the water level element where by virtue of eq. (30) will cause 
the water level to rise. This rise is large for elements close to the 
boundary disturbance and small for elements farther away. This spatial 
change affects the horizontal momentum equation (33) and (34) and induce 
a horizontal flow in the water level elements which in turn tends to re- 
duce the water level through eq. (30). A reverse situation takes place 
when the boundary velocity is decreased. The mass imbalance will result 
in a negative fluxd motion, designated as "downwelling". This will 
cause the lowering of the water level with an attendant horizontal flow 
into the element to increase the water level. 


Boundary Conditions 

The water body boundary surfaces may be classified as follows: 

Type (1): The surface boundary located at the water level free surface. 

Type (2) : The lateral solid-fluid boundary located at the interface be- 
tween the impoundment wall and the water body. 

Type (3): The bottom boundary located at the bottom of the impoundment. 

Type (4); The fluid-fluid boundary located at the interface between the 
water body in the region of interest and surrounding waters. 

For the surface boundary, the pressure is considered atmospheric at the 
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air-water interface and is set equal to zero. For the lateral solid- 
fluid and bottom boundaries, application of the momentum equation in the 
direction normal to the boundaries gives the following equation respec- 
tively 



2 \ 1 _ 


(39) 
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For the fluid-fluid boundaries, the pressure distribution at the inflow 
interfaces is considered to be static. The pressure distribution at any 
outflow interface node is set equal to that of the adjacent node in the 
flow field. The surrounding waters entering a water body is character- 
ized by a uniform channel flow. The pressure variation at any cross- 
section of a uniform channel flow may be proved to be hydrostatic by 
applying the Bernoulli's equation. This equation is applied to two 
streamlines, one at the top of the channel and the other at an arbitrary 
depth. Since, in steady, frictionless, one-dimensional, uniform channel 
flow the constant in the Bernoulli equation is the same for all stream- 
lines, it can be easily shown that the pressure distribution across the 
channel, as well as the interface between the channel and the impound- 
ment is hydrostatic. 


For the surface boundary, the wind effect produces shear at the water sur- 
face which is equal to the shear stress in the fluid near the surface as 
indicated in eqs. (33) and (34) based on the assumptioii that the wind 
velocity components are steady and moderate with zero vertical components 
and the skin coefficients may be taken from experimental measurements 
[18, 13]. For the lateral solid-fluid and bottom boundaries, the velo- 
city components normal to the boundaries is always zero. The tangential 
velocity components is calculated by the application of no-slip condition. 
With this condition, both tangential velocity components at the solid 
boundary are set equal to zero. The no-slip boundary condition is ideal 
for very fine grid spacing. However, computational economy requires a 
coarse grid spacing for most dynamic three-dimensional hydrothermal analy- 
ses. Since, a coarse grid does not afford sufficient resolution at the 
boundary, a limited amount of numerical error will be induced in the solu- 
tion. For the fluid-fluid boundaries, the velocity distribution at the 
inflow interfaces is set equal to the channel inflow velocity distribution 
which is considered to be a known function of time. The velocity dis- 
tribution at any outflow interface node is set equal to that of the 
adjacent node in the flow field. 


For the surface boundary, the rate of heat dispersion through the water 
(involving both conduction and turbulent thermal diffusion) is equal to 
the heat transfer from water to air as indicated in eqs. (32) and (35) 
based on values of the heat exchange coefficient K and equilibrium water 
temperature E,both estimated in terms of meteorological conditions [19]. 


ANN 


- 12 - 



VI-B-49 


For the lateral solid-fluid and bottom boundaries the heat transfer is 
considered to be negligible leading to zero temperature gradient at the 
boundary. The temperature distribution at the fluid— fluid boundaries is 
treated similar to the velocity distribution. The temperature at the 
inflow interface is set equal to the channel inflow temperature which is 
considered to be a known function of time. The temperature at any out- 
flow interface node is set equal to that of the adjacent nodes in the 
flow field. 


The water level at the fluid-fluid boundaries is treated similar to the 
temperature distribution. The water level at the inflow interface is set 
equal to the channel inflow water level which is considered to be known. 
The water level at outflow interface nodes are set equal to that of the 
adjacent nodes in the flow field. For the lateral solid-fluid, the water 
level is calculated from vertical velocity components at surface boundary 
nodes . 


* * * * 

In the development of pressure equation, values of dQ^/dx 9Qy/9y , and 
90*/ 9z* should be evaluated in the flow field. It was shown earlier that 
variables Qy> and are related to the velocity field by eqs. (24), 
(25), and (26). It should be realized that the only spatial derivatives 
of Qx needed in the present analysis is 9Q*/9x*. Calculation of 9Q*/9x* 
in the fluid cell next to the boundary requires the calculation of 
Q* at the boundaries normal to the x axis. Similarly, calculation of 
9Q*/9y* and aQz/Sz* require the calculation of Qy and Qg at the bound- 
aries normal to the y and z axis respectively. For the surface boundary, 
only Q* is needed. As stated earlier, since the vertical velocity com- 
ponent in the water level element is small, eq. (26) gives 

(q”) * * = - 

z = z 

s 0 


(41) 


For the lateral solid-fluid boundaries, only the component of Q normal 
to the solid boundary is required. Since the velocity component normal 
to the solid boundary is zero, eq. (25) will, therefore, give 
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For the bottom boundary only Q is needed. Simplification of eq. (26) 
at the bottom boundary gives 
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(43) 


For the fluid-fluid boundary, the surrounding waters entering the water 
body is characterized by a uniform channel flow, with a uniform velocity 
at any cross section and zero axial pressure gradient. As an example, ap- 
plication of the momentum equation along the channel flowing in the x- 
direction gives 
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(Q*) * = 0 (44) 

^ X = 0 

NUMERICAL SOLUTION 

The numerical scheme, used in this study, is mainly based on the applica- 
tion of spatial and temporal finite difference technique to the governing 
equation described in the previous section. This scheme is divided into 
the following steps [17]: 

1) Initialization and setting of boundary conditions 

2) Calculation of vertical velocity component 

3) Calculation of pressure distribution 

4) Calculation of the time derivatives for the horizontal velocity com- 
ponents and temperature 

5) Calculation of the time derivative for the water level, horizontal 
velocity components and temperature in the water level elements 

6) Time integration and updating of the horizontal velocity components 
and temperature 

7) Time integration and updating of the water level, velocity components, 
and temperature in the water level elements 

8) Calculation of density distribution 

After the completion of the last step, the problem time is incremented, 
the computational procedure is repeated starting at the second step, and 
the integration cycle is continued until the final problem time is 
reached. * 

The initial conditions of the problem consist of the initial values of 
temperature, velocity components, and pressure at time t = 0. These 
quantities are most easily obtained by the application of the present 
dynamic analysis. To demonstrate this application, • let us consider that 
it is desired to determine the dynamic response of a water body when it 
is subjected to a prescribed thermal discharge. The solution to the 
problem can be obtained in two steps. In the first step, the dynamic 
response is obtained by setting the thermal discharge temperature equal 
to the river inflow temperature. Maintaining the atmospheric and the in- 
flow conditions unchanged, the water body will reach equilibriinn con- 
ditions and the pressure, velocity and temperature distribution will be 
obtained. In the second step, these equilibrium conditions are entered 
as the initial conditions and the thermal discharge temperature is set 
equal to the prescribed value. Again, if the atmospheric and inflow 
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conditions are held unchanged, the water body will reach new equilibrium 
conditions. An examination of this dynamic response will enable the pro- 
gram user to predict the effects of the thermal discharge on the state of 
the water body. The results of this study could then be used to deter- 
mine if the state's allowable temperature standards have been violated. 
These standards generally vary from state to state. Many states impose a 
maximum allowable water temperature, a temperature rise, and a rate of 
temperature rise on the use of water bodies. 


From the temporal viewpoint, the boundary conditions stated earlier are 
of two types — time- independent and time-dependent. The setting of the 
time-independent boundary conditions are undertaken prior to the begin- 
ning of* the dynamic solution- The time-dependent boundary conditions are 
satisfied during the execution of the integration cycle. For example, the 
time-dependent boundary conditions on vertical velocity and pressure are 
satisfied during the execution of steps 2 and 3 stated above, respectively. 


Numerical Stability and Convergence 


The numerical procedure, presented in the previous section, involves spat- 
ial integration for the calculation of the vertical velocity components 
and pressure, as well as time integrations for the calculation of the 
horizontal velocity components and temperature. These spatial and time 
integration procedures are plagued by numerical stability problems caused 
by round-off errors in the numerical solution and convergence problems 
caused by the truncation error due to the finite magnitude of the space 
and time increments. A summary of the results for the selection of space 
and time increments are presented hereunder [17]: 


(Ax. + 
1 

(Ay. + 
(Az^ + 




Ax.f-l) ^ 

(45) 


(46) 

4D 


^\-l^ 4 "jw] 

(47) 


At ^1/(8 D^/(Ax^ + + 8Dj^/(Ay^ + Ayj_^)^ 

+ 8Dy(Azj^ + Azj^_^)^) (48) 

At ^l/(2lu|/(Ax^ + ^^i-l^ 2jv|/(Ayj + Ay^^^) ) (49) 


|u|, |v|, and |w| are the estimated maximum absolute values of the velo- 
city components which could occur during the transient. It is true that 
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these values are not known a priori, but an estimate of the maximum values 
can generally be made. 


PRESENTATION OF RESULTS 

The mathematical formulation and the computer program, developed in this 
study, is verified by applying the dynamic analysis to a number of cases 
where the final steady-state solution is either known experimentally or 
the flow pattern has been determined by other investigators. These veri- 
fication studies are considered in the next section. 

Verification Studies 


The following problems were selected for verification purposes: 

1) Experimental studies on the laminar flow development in a square 
duct [20] 

2) Natural circulation in a pond partially heated from the bottom or 
the side 

The first study consists of a duct with a square cross section with a 
spatially and temporally uniform input velocity at the duct entrance. 

Since the geometric dimensions of the laboratory model used in the 
experimental studies were small and our interest here is in water bodies 
of appreciable size, the similarity between the present and the experi- 
mental study is based on two dimensionless parameters; 1) the quotient 
of aspect ratio x/D to Reynolds number; and 2) the- ratio of the fluid 
velocity to the average inlet velocity. The geometric and hydraulic in- 
put data for this case is given in Table I and Fig. 2a. The computational 

grid used in this dynamic analysis consists of 13 x 5 x 5 mesh points 

with non-uniform spacing in the longitudinal direction and uniform spac- 
ing in the square cross section as shown in Table II. The initial 
values of the velocity components in the flow field were selected uniform 
in the longitudinal direction equal to the xnlet velocity with zero 
transversal and vertical velocities. The dynamic solution was continued 
until the flow pattern in the duct reached steady state. 'Typical time 

histories for a number of upstream and downstream points along the duct 

axial central plane are obtained and used to plot the final steady state 
centerline and vertical velocity profiles in the square duct. These 
velocity profiles, in non-dimensional form, are compared in Figs. 2b and 
2c with the experimental results of Goldstein and Kreid [20] who meas- 
ured the laminar flow development in a square duct using laser-doppler 
flowmeter. This comparison shows that the results of this study are in 
good agreement with the experimental data. This confirms that the pre- 
sent formulation can predict the three dimensional flow behavior of water 
bodies . 

The second study consists of natural circulation in a pond heated from the 
bottom or the side with uniform Initial temperature, zero initial velo- 
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cities and zero inlet and outlet mass flux. The geometric and hydraulic 
input data for this case is given in Table III. The computational grid 
used in the dynamic analysis consists of 10 x 6 x 6 mesh points with uni- 
form spacing in the longitudinal and transversal directions and non-uni- 
form spacing in the vertical direction. Two cases were examined. In 
these cases, a step temperature is applied to a portion of side wall or 
the bottom surface of the pond. The dynamic response of the pond is ob- 
tained for each case and the natural circulation flows developed are com- 
pared with expected flow patterns. Typical natural circulation flow 
patterns for partially heated side wall and bottom surface are shown in 
Figs. 3, 4, 5, and 6 at times t = AO and t = 25 seconds respectively into 
the transient. An examination of these flow patterns demonstrates that for 
the case of partially heated side wall one natural circulation vortex in 
transversal direction and two symmetric natural circulation vortices in 
the longitudinal direction are developed which loose their strength as 
the distance from the heated wall increases (see Figs. 3,4). For the case 
of partially heated bottom surface, two symmetric natural vortices in ver- 
tical-longitudinal and vertical-diagonal planes are developed which loose 
their strength with increasing distance from the heated surface (see Figs. 
5,6). In both cases, described above, the natural circulation vortices 
formed and the resultant mixing are caused by density differences in the 
water body. Thesfe results establish that the present formulation can 
predict the three dimensional flow and thermal aspects of water bodies. 

Circulation and Stratification in Water Bodies 


The present mathematical formulation is applied to the study of circula- 
tion and stratification of large water bodies. The water body is con- 
sidered to be flowing initially at a uniform velocity and temperature and 
suddenly exposed to a 90° angle jet with higher velocity but the same tem- 
perature. This problem is referred to as three-dimensional non-buoyant 
jet in a cross current. When the dynamic problem reaches steady state, 
the water jet temperature is suddenly raised to simulate a thermal dis- 
charge. This problem is referred to as three-dimensional buoyant jet in 
a cross current. To facilitate the understanding of the results, the two 
problems indicated above (i.e. non buoyant and buoyant jets) are discussed 
separately but the time histories of the results are shown on the same 
plots for easy reference. 

1) Three Dimensional Non-Buoyant Jet in a Cross Current 

The problem of a three-dimensional non-buoyant jet in a cross current, as 
modeled in Fig. 1, was first analyzed. The geometric and hydraulic input 
data for the case studied are given in Table IV. The water body is as- 
sumed to be flowing initially at a uniform velocity of 0.40 ft/sec when 
suddenly exposed to a jet velocity of 2.00 ft/sec at a 90 degree angle. 

The surface wind velocity is considered to be zero. Furthermore, the 
temperature of the thermal discharge is equal to the water body tempera- 
ture to simulate a non-buoyant jet. A three-dimensional grid was super- 
imposed on the flow field as detailed in Table V. The boundaries of this 
grid coincide with the physical river boundaries. The nodal points, where 
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the hydrothermal variables are defined, are located at point ijk of each 
grid and are shown in Fig. 7. The time histories of the dynamic variables 
are shown, for points A,B,C,D,E.F,G,and H marked in Fig. 7, as follows: 

1) Velocity components and water level for two nodes A and B in front of 

the incoming jet at the free surface shown in Figs. 8 and 9 

2) Velocity components for two nodes C and D in front of the incoming jet 

at 12.75 feet below the free surface shown in Figs. 10 and 11 

3) Velocity components and water level for two nodes E and F located up- 
stream and downstream respectively at the free surface shown in Figs. 12 
and 13 

4) Velocity components of two nodes G and H located upstream and down- 
stream respectively at 12.75 feet below the free surface shown in Figs. 

14 and 15 

Examination of the above plots shows the following features : 

Considering the flow at the water surface, the transversal velocities at 
points A and B, in front of the incoming jet, rise initially with a sub- 
sequent rise in the vertical velocity leading to a rise in water level in 
the neighborhood of the incoming jet which in turn affects the field flow 
conditions as time progresses. Since the total head associated with the 
jet and river flow entries are constant, on a short term basis an increase 
in the water level must be accomplished with a corresponding decrease in 
the field velocity. The disturbance caused by the incoming jet travel 
from the discharge point both upstream and downstream of the river. This 
demonstrates the propagation of the surface gravity waves (surge waves). 

As time progresses, these disturbances reach the river exit and are re- 
flected back upstream. This demonstrates the reflection of the surface 
gravity w^ves from the river exit. However, on a long term basis, a fur- 
ther increase in the water level in the neighborhood of the incoming jet 
causes additional horizontal velocity components for further downstream 
points which will gradually affect the downstream flow. These flow pat- 
terns are observed both in axial and transversal direction as discussed 
hereunder : 

a) The increase in the water level initially decelerates the axial flow 
on a short term basis. With a further increase in the water level, the 
axial flow at various points accelerate to their final steady state val- 
ues which are larger for the downstream points and smaller for the up- 
stream points from the incoming jet position as expected. These effects 
can be clearly observed in u/Ug curves for points E, B and F in Figs. 12, 

9 and 13 where the transient start with a small dip in the flow rate fol- 
lowed by a steep rise settling to its final steady state value. 

b) Similarly, the transversal flow generated by the incoming jet, in- 
creases the water level, which in turn, decelerates the transversal flow. 
However, since, unlike the axial flow patterns, the boundary condition in 
the transversal direction (solid wall) does not accomodate flow, the trans- 
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versa! flow reaches its final steady state value without undergoing large 
swings observed in the axial flow curves. These effects can be clearly ob- 
served in v/Uq curves for points A and B in Figs. 8 and 9 where the trans- 
versal flow after the initial rise undergoes a reduction in magnitude 
followed by a mild rise settling to its final steady state values. 

At the end of transients discussed above, the increase in transversal flow 
caused by the incoming jet is accommodated by an increase in the axial flow. 
For this reason, the water level reaches a steady state value which in turn 
results in zero vertical velocity and constant axial and transversal flows, 
observed in Figs. 8 and 9, as the end of non-buoyant analysis is approached. 

Examining the flow at 12.76 feet below the water surface, the axial flow 
exhibit a pattern similar to that of the water surface except that the re- 
sults are further accentuated due to the bottom boundary condition. Figs. 
11, 14 and 15 show that the flow transients u/Uq start with a relatively 
large dip followed by relatively smaller rise as it approaches its final 
steady state value. The effect of the non-slip bottom boundary condition 
is to reduce the final steady state value to a quantity below that of the 
surface and to enlarge the initial dip. Furthermore, the shear initiated 
by the incoming jet creates transversal flow components v/Uq at this level 
as observed in Figs. 10 and 11. The transversal flow reduces from the top 
to the bottom. At the end of transients discussed above, the increase in 
transversal flow caused by the incoming jet cannot be fully accommodated 
by an increase in the axial flow in cells near the solids boundaries. For 
this reason, a downward velocity component develops which induces a verti- 
cal downward velocity in the upper cells, observed in w/Uq curve in Fig. 

10. This downward flow extends in a few cells from the jet entrance and 
is changed to an upwards flow in cells in the center of the flow field. 

2) Three Dimensional Buoyant Jet in a Cross Current 

The problem of a three-dimensional buoyant jet in a cross current is next 
analyzed. The geometric and hydraulic input data for this case is exactly 
similar to the non-buoyant jet case detailed in Fig. 1 and Table IV. The 
water body is assumed to be initially under the steady state conditions 
reached in the non-buoyant case discussed earlier when the thermal dis- 
charge temperature is suddenly Increased from 75°F to 90°F. 

The dynamic buoyant jet problem is analyzed in a manner similar to the non- 
buoyant jet case. However, in view of the introduction of thermal effects, 
the water body becomes stratified as discussed hereunder. The time his- 
tories of the dynamic variables are shown, for points A,B,C,D,E,F,G, and 
H marked on Fig. 7, as a continuation of the non-buoyant time history plots 
in Figs. 8 to 15. Examination of these plots shows the following features: 

The effect of the heated thermal discharge entering the flow field can be 
decomposed into two components: 1) unheated discharge entering the flow 

field studied in the non-buoyant case problem; and 2) thermal effect of the 
discharge considered in a heated wall studied earlier under the topic of 
ponds partially heated from the side. According to the above decomposition. 
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the velocity time history plots on Figs. 8 to 15 for buoyant case should 
be the sum of the non-buoyant case and the corresponding heated wall prob- 
lem. 

Examining the flow at the water surface at points A and B, in front of the 
incoming jet, the transversal velocity for the buoyant case shows a small 
increase in the steady state as compared with the non-buoyant results due 
to the natural circulation caused by the heated wall. This additional 
transversal velocity induces an axial velocity transient similar in nature 
to the non-buoyant case. However, since the magnitude of the transversal 
velocity transient is small, the resultant axial velocity transient would 
also be small. This behavior can be clearly observed in u/U^ and v/Uq 
curves in Figs. 8 and 9 which also show the temperature transients T/Tq 
for both points A and B. A similar pattern can also be observed at the 
surface points E and F in Figs. 12 and 13 as well as the points C, G, and 
H located at 12.75 feet below the water surface in Figs. 10, 14 and 15. 

At the end of the transients, the coldest and the hottest regions are at 

the river and the thermal discharge entrances respectively. For this rea- 

son, the strongest natural circulation patterns would develop mainly be- 
tween these coldest and hottest regions. A comparison of the surface velo- 
city vectors at points A and E and the velocity vectors at 12.75 feet 
depth for buoyant and non-buoyant cases confirm the existence of this nat- 
ural circulation pattern. The surface velocity field at the end of the 
transient is shown in Fig. 16. The slowing down of the entrance flow near 
the thermal discharge entrance and the turning of the thermal discharge 
and incoming flows as a result of their interaction can be clearly seen 
in this figure. 

The thermal plume effect is shown in Figs. 17 to 22. Isotherms are plot- 
ted for the water surface as well as the vertical planes in the longitudi- 

nal and transversal direction at times 60 and 500 seconds after the ini- 
tiation of the heated discharge. An examination of these plots indicates 
the following features : 

1) The stratification pattern is three dimensional in nature and shows 
the expansion of the thermal plume in all directions. In view of the pre- 
vailing advective currents this expansion is much more intense toward the 
downstream as compared to the other directions (vertical and upstream) as 
expected . 

2) The three distinctive regions which constitute the characteristic be- 
havior of a stratified water body (epilimnlon, thermocline and hypolimnion 
can be clearly observed in Fig. 22. 

3) These results provide the water temperature rise and the rate of tem- 
perature rise needed for the assessment of the extent of thermal pollution 
in the water body. 

As a further verification on the validity of the results a mass and energy 
balance was performed on the water body as shown In Table VI. This table 
shows that conservation equations are satisfied with a good accuracy. The 
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small magnitude of the surface heat transfer verifies the often used as- 
sumption that, at high values of equilibrium temperature used herein, the 
contribution of heat exchange to the atmosphere is insignificant and may 
be neglected in simplified analyses for the purpose of near field studies. 


NOMENCLATURE 


P 

Cfx 

Cfy 

•^o 

Dh 

Dv 

E 


H 

K 

P 

Ro 

So 

To 

T 


t 

Uo 

u 

V 

w 


Wx 

Wy 

x,y,z 

Axi 


4zt 


Specific heat of water 
Skin coefficient along x-axis 
Skin coefficient along y-axis 
Half width of thermal discharge 
Horizontal eddy diffusivity of heat 
Vertical eddy diffusivity of heat 
Water equilibrium temperature 
Froud number 

Gravitional acceleration 

Water level height in the water level element 

Heat exchange coefficient for water-air interface 

Pressure 

Reynolds number 

Stanton number 

Thermal discharge temperature 
Local water temperature 
Time 

Thermal discharge velocity 

Horizontal velocity components in x-direction 
Horizontal velocity components in y-direction 
Vertical velocity components in z-direction 
Wind velocity components in x-directlon 
Wind velocity components in y-direction 
Cartesian coordinate system 
Dimension of element i in x-direction 
Dimension of element j in y-direction 
Dimension of element k in z-direction 


Greek Symbols 


P 

PQ 

Pa 

^h 

Vv 

''O 

2 

V 


Local density of water 
Thermal discharge density 
Air density 

Horizontal eddy viscosity 

Vertical eddy viscosity 

Thermal discharge kinematic viscosity 

2 2 2 

3 ^ 9 ^ 9 ^ 

Laplacian operator — ^ ^ + — r 

3x 3y 9z 


Superscripts 

* Refers to non-dimensional quantities 


-21- 


ANN 



VI-B-58 


Subscripts 

0 Refers to thermal discharge 

s Refers to values of variables at the free surface 

i,j,k Indices referring to the location of nodal points along x,y, 
and z axes respectively 
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TABLE I - THE GEOMETRIC AND HYDRAULIC INPUT DATA 
FOR LAMINAR FLOW IN A SQUARE DUCT 


Specifications 


Dimensions 


British Unit SI Unit 


Water Body Length 

280 ft 

85.34 

m 

Water Body Width 

72 ft 

21.946 

m 

Water Body Depth 

72 ft 

21.946 

m 

Inlet Water Velocity 

1 ft/sec 

0.3048 m/; 

Water Temperature 

75 °F 

23.89 

oc 

Reynolds Number, Rq 
H ydraulic Diameter of 

20.83 

20.83 


Square Duct, D 

72 ft 

21.946 

m 


TABLE II - GRID DIMENSIONS FOR THE LAMINAR 
FLOW IN A SQUARE DUCT 

Element 

No. 


X 2. ■ Z. 



ft 

m 


m 

ft 

m 

1 

10 

3.048 

18 

5.486 

18 

5.486 

2 

10 

3.048 

18 

5.486 

18 

5.486 

3 

10 

3.048 

18 

5.486 

18 

5.486 

4 

10 

3.048 

18 

5.486 

18 

5.486 

5 

20 

6.096 

18 

5.486 

18 

5.486 

6 

30 

9.144 





7 

30 

9.144 





8 

30 

9.144 





9 

30 

9.144 





10 

30 

9.144 





11 

30 

9.144 





12 

30 

9.144 





13 

30 

9.144 
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TABLE III - GEOMETRIC AND HYDRAULIC INPUT DATA 
FOR PARTIALLY HEATED POND 


Specifications 


Dimensions 



British 

Unit 

SI Unit 


Water Body Length 

540 

ft 

164.592 

m 

Water Body Width 

250 

ft 

76.2 

m 

Water Body Depth 

20.5 

ft 

6.248 

m 

Water Body Temperature 

75 

°F 

23.89 

oc 

Water Body Equilibrium 





Temperature 

74 

op 

23.33 

°C 

Temperature of Heated Area 

100 

Op 

37.77 

°C 

Reference Velocity, Ug 

1 

ft/sec 

0.3048 

m/ sec 

Reference Length, dg 

50 

ft 

15.24 

m 

Reference Temperature, Tg 

100 

Op 

37.77 

°C 

Reference Time, tg 

50 

sec 

50 

sec _ 

Reference Density, Pg 

61.9963 

Ibm ft" 

960.590 

kg m 

Reference Pressure, pg 

21.52 

Ibf in~^ 

14639.4 

kg m~2 

Reynolds Number, Rg 

0.6061 

X 10^ 

0.6061 X 

107 

Heat Exchange Coeff., K 

100 

Btu/f t^day°F 

23.64 

W/m2oc 


TABLE IV - GEOMETRIC AND HYDRAULIC INPUT DATA FOR THREE-DIMENSIONAL, 
NON-BUOYANT AND BUOYANT JETS IN A CROSS CURRENT 


Specifications 


Water Body Length 
Water Body Width 
Water Body Depth 
Jet Width, 2 dg 
Jet Depth 
Jet Velocity, Ug 
River Velocity 
Water Body Temperature 
Water Body Equilibrium 
Temperature 
Thermal Discharge 
Temperature, Tg 
Reference Time, tg 
Reference Density, Pg 
Reference Pressure, pg 
Reynolds Number, Rg 
Heat Exchange Coef f . , K 


Dimensions 


British Unit 


600 

ft 

182.88 

m 

360 

ft 

109.73 

m 

16.5 

ft 

5.03 

m 

100 

ft 

30.48 

m 

/ 

5.75 

ft 

1.75 

m 

2 

ft/sec 

.61 

m/sec 

.4 

ft/sec 

.12 

m/sec 

75 

Op 

23.89 

°C 

74 

Op 

23.33 

°C 

90 

Op 

32.22 

°C 

25 

sec 

25 

sec 

62.1156 

Ibm ft~'' 

962.730 

kg m"^ 

21.56 

Ibf in-2 

14672.0 

kg m"2 

0.1212 

X 108 

0.1212 

X 10® 

100 

Btu/ft^day°F 

23.64 

W/m2oc 
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TABLE V - GRID DIMENSIONS FOR THREE DIMENSIONAL 
NON-BUOYANT AND BUOYANT JETS IN A CROSS CURRENT 


Element 


No. 




X 




z 


ft 

m 


m 

it 

m 

1 

65 

19.81 

50 

15.24 

3.5 

1.07 

2 

60 

18.29 

50 

15.24 

4.0 

1.22 

3 

55 

16.76 

55 

16.76 

5.0 

1.52 

4 

50 

15.24 

60 

18.29 
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TABLE VI - 

MASS AND ENERGY BALANCE 



Mass or Energy 
Balance 

Incoming 

Flow 

Thermal 

Discharge 

Outgoing 

Flow 

Surface Percent 
Flow Error 

Mass Balance, 
in 10^ Ibm/sec 

0.12252 

0.08464 

0.20748 

0.0 

0.15 

Energy Balance, 
in 10° Btu/sec 

0.09189 

0.07454 

0.16204 

0.15096 

X 10“^ 

0.18 
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FIG I. SCHEMATIC DIAGRAM OF THE WATER BODY 
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(a) THE SQUARE DUCT CONFIGURATION 



(b) CENTER-LINE VELOCITY DEVELOPMENT 



0 2 

(c) DEVELOPMENT OF VELOCITY PROFILE 
CENTRAL PLANE 


FIG 2. VELOCITY DEVELOPMENT IN A SQUARE DUCT 
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O’URAL CIRCULATION AT y = IOOFT. AND t=40 SEC 
A POND PARTIALIY HEATED FROM SIDE. 



HEATED AREA 
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FIG 4. NATURAL CIRCULATION AT X* 240 FT. 
AND t = 40 SEC. IN A POND PARTIALLY 
HEATED FROM SIDE. 
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FIG 5. NATURAL CIRCULATION AT y = IOO FT. AND t= 25 SEC. 
IN A POND PARTIALLY HEATED FROM BOTTOM. 



20.5 FT 




FIG 6. NATURAL CIRCULATION AT t=25SEC IN VERTICAL 
DIAGONAL PLANE IN A POND PARTIALLY HEATED 
FROM BOTTOM 
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FIG 7. GRID WORK WITH VARIABLE MESH SIZE 
SUPERIMPOSED ON THE WATER BODY 
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DIMENSIONLESS VARIABLES 
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FIG 8. TIME HISTORIES OF VARIABLES AT POINT "A" IN THE WATER BODY 
SUBJECTED TO NON-BUOYANT AND BUOYANT JETS 
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DIMENSIONLESS VARIABLES 



FIG 9. TIME HISTORIES OFVARIABLES AT POINT "B“ IN THE WATER BODY 
SUBJECTED TO NON-BUOYANT AND BUOYANT JETS 
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DIMENSIONLESS VARIABLES 



FIG 10. TIME HISTORIES OF VARIABLES AT POINT "C" IN THE WATER BODY 
SUBJECTED TO NON- BUOYANT AND BUOYANT JETS 
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DIMENSIONLESS VARIABLES 
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FIG II. TIME HISTORIES OF VARIABLES AT POINT "D" IN THE WATER BODY 
SUBJECTED TO NON -BUOYANT AND BUOYANT JETS 
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DIMENSIONLESS VARIABLES 



FIG 12. TIME HISTORIES OF VARIABLES AT POINT "E" IN THE WATER BODY 
SUBJECTED TO NON -BUOYANT AND BUOYANT JETS 
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DIMENSIONLESS VARIABLES 



FIG 13. TIME HISTORIES OF VARIABLES AT POINT “f“ IN THE WATER BODY 
SUBJECTED TO NON-BUOYANT AND BUOYANT JETS 
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FIG 14. TIME HISTORIES OF VARIABLES AT POINT "G" IN THE WATER BODY 
SUBJECTED TO NON-BUOYANT AND BUOYANT JETS 
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DIMENSIONLESS VARIABLES 
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FIG 15. TIME HISTORIES OF VARIABLES AT POINT "H" IN THE WATER BODY 
SUBJECTED TO NON-BUOYANT AND BUOYANT JETS 
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FIG 16. SURFACE VELOCITY FIELD AT t= 500 SEC. IN THE WATER 
SUBJECTED TO A BUOYANT JET 
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FIG 18. VERTICAL ISOTHERMS AT y = 100 FT. AND t = 60 SEC. IN THE ^ 
WATER BODY SUBJECTED TO A BUOYANT JET 
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FIG 20. SURFACE ISOTHERMS AT t=500 SEC. IN THE WATER BODY 

SUBJECTED TO A BUOYANT JET | 
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AND SOO SEC. IN THE WATER BODY 
SUBJECTED 10 A BUOYANT JET 
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